The second Hopf bifurcation in lid-driven square cavity
Wang Tao1, 2, Liu Tiegang1, †, Wang Zheng3
LIMB and School of Mathematical Science, Beihang University, Beijing 100191, China
School of Mathematics and Information Science, North Minzu University, Yinchuan 750021, China
Wuhan Maritime Communication Research Institute, Wuhan 430079, China

 

† Corresponding author. E-mail: liutg@buaa.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11601013 and 91530325).

Abstract

To date, there are very few studies on the second Hopf bifurcation in a driven square cavity, although there are intensive investigations focused on the first Hopf bifurcation in literature, due to the difficulties of theoretical analyses and numerical simulations. In this paper, we study the characteristics of the second Hopf bifurcation in a driven square cavity by applying a consistent fourth-order compact finite difference scheme recently developed by us. We numerically identify the critical Reynolds number of the second Hopf bifurcation located in the interval of (11093.75,11094.3604) by bisection. In addition, we find that there are two dominant frequencies in its spectral diagram when the flow is in the status of the second Hopf bifurcation, while only one dominant frequency is identified if the flow is in the first Hopf bifurcation via the Fourier analysis. More interestingly, the flow phase portrait of velocity components is found to make transition from a regular elliptical closed form for the first Hopf bifurcation to a non-elliptical closed form with self-intersection for the second Hopf bifurcation. Such characteristics disclose flow in a quasi-periodic state when the second Hopf bifurcation occurs.

PACS: ;05.45.Pq;;05.70.Jk;
1. Introduction

The two-dimensional (2D) flows in driven cavities have been extensively studied in the past 60 years,[17] starting from the pioneering work by Kawaguti. Among them, several works have contributed to establish a clear picture of the steady solutions to the problems for Reynolds numbers up to several thousands, mainly by employing vorticity–velocity and biharmonis solvers.[2,4,6]

In recent decades, scientists have shifted their research focus on studying the first Hopf bifurcation,[810] in which the flow transfers its stationary state to non-stationary periodic state. In literature, there are majorally two kinds of methods applied to study the first Hopf bifurcation. One is through theoretical analyses, the other is via numerical simulations. The first attempt was made by Hou et al.,[8] who numerically predicted that the critical Reynolds number of Rec1 was larger than 7500 for the occurrence of the first Hopf bifurcation with a lattice Boltzmann method. Immediately, Poliashenko and Aidun[9] confirmed the prediction of Hou et al. and gave out Rec1 = 7763 ± 2% by a low order finite element method. It seems true that the critical Reynolds number of the first Hopf bifurcation is beyond 7500 both from theoretical analysis and numerical simulation. In theoretical analysis, Fortin et al.[11] predicted the critical Reynolds number of Rec1 = 8000 for the first Hopf bifurcation, through an eigenvalue analysis of the linearized Navier–Stokes equations by means of a finite element spatial discretization. Late, Cazemier et al.[12] found Rec1 was at 7819.0 through the eigenvalue analysis of the linearized Navier–Stokes equations but with a proper orthogonal decomposition to solve the linearized system. On the numerical side. Auteri et al.[13] applied a singularity subtraction technique and second-order spectral projection method to locate the Rec1 was in the interval (8017.6,8018.8). They claimed that an asymptotic solution existed when Reynolds number was less than Rec1 and the asymptotic solution would become periodic when Reynolds number was larger than Rec1. More recently, Bruneau et al.[10] calculated Rec1 was in the interval (8000,8050) by a third-order Murman-like scheme.

With the increase of Reynolds number, the flow experiences transition from the first Hopf bifurcation to the second Hopf bifurcation, where the flow is supposed to be under a quasi-periodic state. Due to the second bifurcation is a further bifurcation of flow instability, it is very difficult to get useful information under the frame of eigenvalue analysis of the linearized Navier–Stokes equations. To date, one can only use numerical simulations to make investigation and very few results disclosed. Auteri et al.[13] numerically forecasted the second Hopf bifurcation with a critical Reynolds number, Rec2, located in the interval of (9685,9765). Because the flow instability starts from the wall boundary layer, where it becomes very thinner under a high Reynolds number. To predict the second bifurcation, the boundary layer has to be resolved accurately with either very fine mesh or high order numerical method. This makes the numerical simulation very costly and challenging. In even worse case, recently, it was found by us[14] that a high order method might suffer numerical instability inside the wall boundary layer if a lower order treatment of boundary conditions is equipped. Such lower order techniques of wall treatment are very popular in applications of high order numerical methods. To fix this problem, we have developed a consistent fourth-order compact finite difference scheme,[14] in which both the order of accuracy and major truncation error term are kept the same for the inner and boundary schemes. Such a high order scheme can efficiently suppress the numerical instability insider the wall boundary layer.

In this work, we will apply the consistent fourth-order compact finite difference scheme to study the characteristics of the second Hopf bifurcation in the driven square cavity. In addition, we employ bisection to locate the critical Reynolds number and analyze the characteristics of the second Hopf bifurcation with Fourier analysis and phase diagram. We discover that the flow phase portrait makes transition from a regular elliptical closed form for the first Hopf bifurcation to a non-elliptical closed form with self-intersection when the second Hopf bifurcation occurs, indicating a quasi-periodic flow behavior.

2. Model description and numerical methods
2.1. Governing equations

In this paper, the 2D lid–driven cavity problems are considered and governed with the 2D incompressible Navier–Stokes equations given in the non-dimensional vorticity–stream function in the form as

where u and v are the velocity components in x and y directions, respectively, Re is the Reynolds number and u = ψy, v = − ψx. Equation (1) is referred to as the stream function equation, equation (2) is referred to as the vorticity equation.

The computational domain is set to {(x,y)|0≤ x ≤ 1,0≤ y ≤ 1}. The boundary conditions are given as follows:

Because the vorticity on the wall is unknown, in the process of solving the Eq. (2), the condition of vorticity on the boundary must be supplemented. In this paper, the numerical boundary conditions of vorticity in Ref. [15] are used, namely,

where the subscripts 1, 2, 3 represent the boundary points, its neighbor, and sub-neighbor points, respectively. Vw is the tangential velocity of the wall. On the sliding wall, Vw = 1. On the solid wall, Vw = 0.

In Fig. 1, the typical flow pattern in a driven cavity is illustrated, where there is a primary vortex (PV) at the center, three corner vortexes locate respectively at the bottom (bottom left vortex (BL), bottom right vortex (BR), and the top (top left vortex (TL)).

Fig. 1. The lid-driven cavity problem: (a) boundary condition and (b) the schematic diagram.
2.2. Numerical methods

The numerical method employed in this work is the consistent fourth-order compact finite difference scheme developed by Wang et al.[14] One can refer to Ref. [14] for details. After the spatial discretization with the consistent fourth-order compact finite difference scheme, an explicit third-order TVD Runge–Kutta scheme (for details see Ref. [16]) was used to do time discretization.

Mesh convergency for the numerical method has been carried out with a driven square cavity with Re = 10000. Table 1 shows the amplitude values and periodic lengths of the u-velocity when the fluid reaches the periodic state (t = 2000) using three uniform meshes (65 × 65,129 × 129 and 257 × 257) at Re = 10000. It can be seen that the results under grid 129 × 129 are in good agreement with those under fine grid 257 × 257 and the maximum error is less than 1%.

Table 1.

Comparison of the amplitude values and periodic lengths of the u-velocity when the fluid reaches the periodic state at Re = 10000.

.

Based on the mesh convergences test of the above benchmark simulation, in this paper, we take the time step size of Δt = 1 × 10−3 and solve the system with 129 × 129 uniform grids in the numerical simulations.

3. Results and discussion

We determine the different states of the flow according to the time evolution of velocity at a fixed point (x,y) = (0.5,0.5) of the cavity.

3.1. The critical Reynolds number of the second Hopf bifurcation

From the introduction, it is known that the flow in a driven square cavity becomes periodic when the Reynolds number is larger than a critical number of Rec1, which is beyond 7500. We confirmed that the flow at Re = 10000 is still periodic in its first Hopf bifurcation as done by Tian et al.[17] We calculated the flow with a higher Re = 12500 for t from 0 to 3000 as shown in Figs. 2(a) and 3(a) for the evolution of u-velocity and v-velocity at the geometric center point. We give the enlarged Figs. 2(a) and 3(a) for t from 2960 to 3000 and found that the flow was under quasi-period as shown in Figs. 2(b) and 3(b) for the evolution of u-velocity and v-velocity at the geometric center point. From the enlarged figures, one can see that the flow from the initial unstable state converges eventually into a quasi-periodic state.

Fig. 2. The time history of u-velocity for Re = 12500: (a) total time; (b) local time.
Fig. 3. The time history of v-velocity for Re = 12500: (a) total time; (b) local time.

As a result, the critical Reynolds number of Rec2 is located in between Re = 10000 and Re = 12500. Now, we find the Rec2 via bisection. Starting from a0 = 10000, b0 = 12500, if at Re = (ak + bk)/2, we get periodic solution, then we let ak=(ak + bk)/2; If at Re = (ak + bk)/2, we get quasi-periodic solution, then we let bk=(ak + bk)/2. We repeat this process and carry out the dichotomy for 12 times, and then we obtain a very small interval. Finally, in this way, the Rec2 has been localized in the small interval (11093.75, 11094.3604). Table 2 shows the process of localizing the range of the Rec2 by bisection.

Table 2.

The process of bisection for confirming Rec2.

.

In a series of Figs. 4 and 5 we present the time history of u(t) for different Reynolds numbers. In these graphes, we can clearly see that during the process of determining the critical Reynolds number, the function u(t) solution is either periodic or quasi-periodic.

Fig. 4. The time history of u-velocity for different Reynolds numbers: (a1)–(f1) total time, (a2)–(f2) local time.
Fig. 5. The time history of u-velocity for different Reynolds numbers: (a1)–(f1) total time, (a2)–(f2) local time.
3.2. The analysis of second Hopf bifurcation

We analyze the flow behavior from the time history of flow velocity, its spectrum and phase diagram for two specific Reynolds numbers of Re = 10000, 12500. For the former, the flow is in the status of its first Hopf bifurcation, for the latter the flow is in its second Hopf bifurcation. Since the properties of u-velocity and v-velocity are the same, we only look at the behavior of the u-velocity. In this way, we may have a clear understanding of the nature of the second bifurcation.

Figures 6(a) and 6(b) show their time histories of u-velocity for the respective Re = 10000 and Re = 12500; the solutions clearly exhibit periodic behaviors with non-dimensional periodic lengths of τ = 1.818 for the former and τ = 1.429 for the latter.

Fig. 6. The local time history of u-velocity for (a) Re = 10000 with τ = 1.818 and (b) Re = 12500, τ = 1.429.

The flow behaviors looks very different from Figs. 6(a) and 6(b). We took Fourier analysis from t = 2960 to t = 2970. The results are shown in Figs. 7(a) and 7(b) for the power spectrum density of u(t) for the respective Re = 10000 and Re = 12500. There is only one frequency of f = 0.55005 found for the flow at Re = 10000, while there are two distinct frequencies for the flow at Re = 12500, one with a frequency of f1 = 0.38654 and the other with a frequency of f2 = 0.77309. This show that the flow at Re = 10000 is periodic, while it is in the quasi-periodic behavior at Re = 12500.

Fig. 7. The u-velocity power spectrum density of u(t) for (a) Re = 10000 and (b) Re = 12500.

We further plotted the phase portraits by using the delay coordinates u(t) and u(t + τ) with τ = 1.818 for Re = 10000 and τ = 1.429 for Re = 12500, respectively. The results are shown in Figs. 8(a) and 8(b). One can observe that a perfect elliptic cycle formed for Re = 10000, while a non-ellipse closed pattern formed with self-intersection for Re = 12500. These features further confirms that the flow is in stable period under the first Hopf bifurcation and in quasi-period for the second Hopf bifurcation.

Fig. 8. Phase portrait of u(t) for (a) Re = 10000, τ = 1.818 and (b) Re = 12500, τ = 1.429.

In all, when Re ∈ (10000,Rec2), the flow state is periodic, which is similar to that for Re = 10000 although the periodic value is different for different Re. And when Re ∈ (Rec2,12500), the flow state is quasi-periodic, which is similar to that for Re = 12500 although the two frequencies are different for different Re.

Finally, we plot the findings in Fig. 9 to give a qualitative summary. The flow is stationary below Rec1, periodic in Re ∈ (Rec1, Rec2), and quasi-periodic in Re ∈ (Rec2,12500).

Fig. 9. Bifurcation diagram.
4. Conclusion

In this paper, we investigated the flow properties under its second Hopf bifurcation in a driven square cavity by applying a consistent fourth-order compact finite difference scheme. We have located the critical Reynolds number of the second Hopf bifurcation in a very narrow interval of (11093.75,11094.3604). Moreover, by analyzing the spectrum and phase portrait, we found that there are two dominant frequencies and the phase portrait changes from an ellipse closed form to the non-ellipse closed form, when the flow is in the status of the second Hopf bifurcation. We conjecture that the flow might experience further bifurcation before it becomes turbulent. This will be our future work.

Reference
[1] Kawaguti M 1961 J. Phys. Soc. Jpn. 16 2307
[2] Erturk E 2006 Int. J. Numer. Meth. Fluids 50 421
[3] Das M K Kanna P 2007 International Journal of Numerical Methods for Heat and Fluid Flow 17 799
[4] Erturk E 2009 Int. J. Numer. Meth. Fluids 60 275
[5] Tang J S Zhao M H Han F Zhang L 2011 Chin. Phys. 20 020504
[6] Nuriev A N Egorov A G Zaitseva O N 2016 Fluid Dyn. Res. 48 61405
[7] Fang L Ge M W 2017 Chin. Phys. Lett. 34 030501
[8] Hou S Zou Q Chen S Doolen G Cogley A C 1995 J. Comput. Phys. 118 329
[9] Poliashenko M Aidun A 1995 J. Comput. Phys. 121 246
[10] Bruneau C H Saad M 2006 Comput. Fluids 35 326
[11] Fortin A Jardak M Gervais J Pierre R 1997 Int. J. Numer. Method Fluids 24 1185
[12] Cazemier W Verstappen R Veldman A 1998 Phys. Fluids 10 1685
[13] Auteri F Parolini N Quartapelle L 2002 J. Comput. Phys. 183 1
[14] Wang T Liu T G 2019 Numer. Math.: Theory, Methods and Applications 12 312
[15] Spotz W F 1998 Int. J. Numer. Method Fluids 28 737
[16] Shu C W Osher S 1988 J. Comput. Phys. 77 439
[17] Tian Z F Liang X Yu P 2011 Int. J. Numer. Method Eng. 88 511